The Ising model in a Bak-Tang-Wiesenfeld sandpile 
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We study the spin-1 Ising model with non-local constraints imposed by the Bak-Tang-Wiesenfeld 
sandpile model of self-organized criticality (SOC). The model is constructed as if the sandpile is 
being built on a (honeycomb) lattice with Ising interactions. In this way we combine two models that 
exhibit power-law decay of correlation functions characterized by different exponents. We discuss 
the model properties through an order parameter and the mean energy per node, as well as the 
temperature dependence of their fourth-order Binder cumulants. We find (i) a thermodynamic phase 
transition at a finite Tc between paramagnetic and antiferromagnetic phases, and (ii) that above Tc 
the correlation functions decay in a way typical of SOC. The usual thermodynamic criticality of the 
two-dimensional Ising model is not affected by SOC constraints (the specific heat critical exponent 
a « 0), nor are SOC-induced correlations affected by the interactions of the Ising model. Even 
though the constraints imposed by the SOC model induce long-range correlations, as if at standard 
(thermodynamic) criticality, these SOC-induced correlations have no impact on the thermodynamic 
functions. 

PACS numbers: 05.65.-|-b; 45.70.Ht; 75.10.Hk; 



I. INTRODUCTION 

The phenomenon of self-organized criticality (SOC) 
attracts a lot of interest in various branches of science 
(see 0, for reviews). Its most intriguing feature re- 
sides in the "spontaneous" formation of scale-free spatio- 
temporal patterns with power-law correlations between 
various quantities. On the one hand, these correlations 
closely resemble those that appear at critical points in 
continuous phase transitions. On the other hand, while 
the critical state in phase transitions is temperature- 
or external-field-driven, a SOC state is thought to form 
without fine-tuning of any external parameters. The re- 
lation between the "classical" and "self-organized" criti- 
cality has been studied by many authors 0,13,0) 01 with a 
conclusion that the origin of the self-organized criticality 
is a continuous absorbing-state phase transition to which 
a SOC system is attracted. 

Perhaps the best-known and most extensively studied 
model of thermodynamic criticality is the Ising model of 
ferromagnetism. Of similar importance for self-organized 
criticality is the Bak-Tang-Wiesenfeld (BTW) model 
|3| of a sandpile. Owing to its elegant mathematical 
structure and many rigorous results 0j H, IS the 
BTW model has practically become a paradigm for self- 
organized criticality studies. 

Specific SOC features of the BTW model, which will 
attract our attention, include the fact that in the BTW 
model correlation functions decay with the distance r as 
r~2d^ where d is the space dimensionality [TllIT^. which 
resembles, but differs from the r~*^'^~^+''^ decay in stan- 
dard critical phenomena. As the behavior of a correlation 
function is a key ingredient of criticality, we will use this 
property to build and analyze a model which, by con- 



struction, is expected to exhibit long-range correlations 
induced by both classical and self-organized critical phe- 
nomena. In defining such a model we follow the idea of 
Ref. 13^: take a standard model of statistical physics 
and significantly reduce its phase space to that of a cor- 
responding SOC model. 

One of the main results of ref. 'l3l, where a combi- 
nation of the Potts and BTW models was investigated, 
is that in the limit of a vanishingly small temperature 
their hybrid system looses its 'self-organized criticality', 
i.e. the power-law correlations disappear. Thus it seems 
worthwile to investigate a hybrid model, numerically, for 
a wide range of temperatures. In particular, we want to 
check if the self-organized criticality imposed by nonlo- 
cal constraints can be detected in the behavior of ther- 
modynamic functions. However, as it was pointed out in 
Ref. ^3 ) because of the non-local nature of the applied 
constraints, analysis of this type of models is a difficult 
problem, still at its infancy, so it is of interest to check 
whether standard concepts and theorems of statistical 
physics, like existence of the thermodynamic limit, er- 
godicity, fluctuation-dissipation theorem or the central 
limit theorem apply to such hybrid systems. 

In our study we choose the two-dimensional spin-1 
Ising model as the Hamiltonian system and the BTW 
sandpile model as the SOC component. The resulting 
model is an equilibrium model that takes (short-range) 
interactions from the Ising model and the (non-local) 
constraints from the BTW model. Note, however, that 
while the Ising model focuses on spin configuration at 
thermodynamical equilibrium, the BTW model attempts 
to describe a nonequilibrium, dynamical process where 
grains of sand are continually added to the lattice to form 
"avalanches" of different sizes and durations. 

The construction of such a "hybrid" model employs the 
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fact that both component models are built on a lattice 
with special variables ("spins" Sj or "grain heights" hj) 
attached to each lattice node, and the values of these 
variables are limited to a few integers, e.g. sj G {—1, 0, 1} 
or hj g {0,1,2}. This property allows us to map the 
states of one model onto the states of the other one. To 
reduce to the minimum the number of possible mappings 
of the spins sj of the Ising model to the sand heights hj 
of the BTW sandpile, we choose a two-dimensional (2D) 
honeycomb lattice (coordination number z = 3). 

In a hybrid system like this the high-temperature limit 
is determined entirely by the constraints of the SOC 
component; consequently, the spin correlation function 
should decay with the distance r as r~'^ [Til IT^ . On the 
other hand, the ground state is expected to be controlled 
mainly by the Hamiltonian, and at sufficiently low tem- 
peratures interactions should force the system into an 
ordered phase. Consequently, one can expect a thermo- 
dynamic phase transition at some finite temperature Tc. 
For the standard 2D Ising model without an external field 
this transition is continuous and the spatial correlations 
decay as with rj — 1/4 JJJ. Therefore, if introduc- 
tion of the SOC constraints turns out to be too weak 
to change the nature of this transition, the hybrid model 
should exhibit a power-law decay of correlation functions 
both for T — !■ cxD and at Tc. 

Although the Ising model has been modified in many, 
many ways, most of the modifications are of local char- 
acter. These include building the model on fractal 
or small-world llal lattices as well as introduction of 
quenched 0, Il8|| or kinetic 0, |23| disorder. In our 
present approach the modification is nonlocal - one has 
to scan the whole lattice to decide whether a given con- 
figuration is allowed or not. This introduces several com- 
plications, but at the same time opens new directions of 
research. In our analysis we shall consider both the ef- 
fects of introducing interactions and temperature to a 
SOC model and the changes introduced by a SOC-like 
constraints into a Hamiltonian system. 

The paper is organized as follows. In section^ we de- 
fine the model. Then, in section UTTl we discuss problems 
that appear in Monte Carlo simulations and we present 
some methods that we have applied to overcome them. 
The results are presented in section HVI They include a 
numerical analysis of the finite-size effects imposed by the 
SOC subsystem, an analytical study of the influence of 
the sandpile model on the number of phases in the Hamil- 
tonian subsystem, and a numerical analysis of the impact 
of the SOC and Hamiltonian components on the critical 
properties of the system. It also contains results of sev- 
eral tests carried out to verify if the hybrid model with 
nonlocal constraints can be treated with standard meth- 
ods of statistical physics. These include tests of whether 
fluctuations of the internal energy can be used to deter- 
mine the speciflc heat and whether far from the critical 
point the central limit theorem can be applied to pre- 
dict the distribution of the energy fluctuations. Finally, 
section is devoted to conclusions. 




FIG. 1: An example of an allowed BTW configuration on a 
honeycomb lattice of linear size L = 3. Empty, shadowed, 
and filled circles represent the nodes with heights hj = 0, 1, 2 
(or spins Sj = 0,-1,-1-1), respectively. The arrows show an 
example of a pair of boundary nodes that interact via the 
Ising Hamiltonian Q. 

II. MODEL 

We define the model on a two-dimensional honeycomb 
lattice of linear size L. An example of such a lattice with 
L = 3 is shown in Fig. ^ The lattice has a shape of a 
big hexagon made of 3L(L — 1) -I- 1 small hexagons and 
it has N = nodes. Most of them are interior nodes 
connected to 3 nearest neighbors, but 6L nodes lie on the 
edge of the lattice and have only 2 connections. 

Each node j can be in one of three states, which can 
be interpreted either as a spin variable Sj G {— 1, 0, 1} or 
the number of sand grains hj G {0, 1, 2}. Of 6 mappings 
of heights hj on spins Sj, we choose the one in which 
hj = 0,1,2 corresponds to Sj — 0,-1,-1-1, respectively 
(actually, due to the Z2 symmetry, the number of non- 
equivalent mappings reduces to 3). 

The Hamiltonian of the model is simply that of the 
spin-1 Ising model, 

H = -jJ2s.Sj (1) 
(hi) 

where J is a real parameter, Sj = 0, ±1, and the sum is 
to be taken over all pairs of nearest neighbors. To mini- 
mize finite-size effects, we have adopted periodic bound- 
ary conditions. To this end we assume that, as is illus- 
trated in Fig. n each boundary node interacts, via iQ), 
not only with its two direct lattice neighbors, but also 
with the node at the location symmetric about the lat- 
tice center. This choice ensures that the ground state in 
the antiferromagnetic case ( J < 0) is doubly degenerated 
and, together with the several consecutive excited states, 
obey the Z2 symmetry. 

In the standard spin-1 Ising model the phase space 
is defined by all possible spin configurations. To each 
such configuration 77 one can assign a non-vanishing 
temperature-dependent probability oc exp[— _ff (77)/A:bT], 
and their number is exactly 3^. However, the states 
of the BTW model, which is a dynamic, nonequilibrium 
model of a sand pile in a continual flux of sand grains, 
have completely different properties. It turns out that all 
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configurations can be divided into two categories: disal- 
lowed and allowed. The former can never be found in the 
steady state^while the latter can be found with the same 
probability 0,11 0. Moreover, the number of allowed 
configurations is with a < 3 8]; hence the number of 
allowed configurations is only a fraction of all configura- 
tions and this fraction tends to as L ^ cxd. 

The main idea behind building a hybrid Hamiltonian- 
SOC model is to take a standard model of equilibrium 
statistical physics, e.g. the spin-1 Ising model, and sig- 
nificantly reduce its phase space to that of a correspond- 
ing SOC model, e.g. the BTW model. Solving a hybrid 
model is then equivalent to finding the partition function 

Z = Y,^iv)eM-H{v)/kBn (2) 
n 

where the sum runs over all states rj of the Hamiltonian 
component, 8 is the characteristic hmction of a SOC 
state, 

9(?7) = 1, if is allowed 
9(77) = 0, if 77 is disallowed 

and fcs is the Boltzmann constant. 

To complete the definition of the model, it remains to 
explain how to distinguish an allowed configuration from 
a disallowed one in the BTW model. A complete dis- 
cussion of this problem is given by Dhar in j3, Q ; here 
we will give only a brief summary of his results: (i) a 
configuration with all hj = 2 is allowed, (ii) If we take 
an allowed configuration {hj} and increase the value of 
hj at some j by one, this hj may exceed the maximum 
value 2; however, if we then relax such an "unstable" 
state (in the way defined below), we will always arrive 
at an allowed configuration, (iii) All allowed configura- 
tions can be reached by iteratively applying step (ii) to 
any allowed configuration. Relaxation of unstable states 
is carried out as follows: for any site j at which hj > 2, 
replace hj with hj — 3 and for every neighbor k of j 
increase hk by 1. Relaxation of any state can be car- 
ried out in arbitrary order, consists of a finite number of 
steps that form the so called avalanches and the resulting 
stable state is unique 0, Q . An example of an allowed 
configuration is shown in Fig. ^ Note that while we use 
periodic boundary conditions for the Hamiltonian part of 
the model, we employ open BC for the relaxation process 
so that excessive "sand grains" can leave the system. 

III. NUMERICAL IMPLEMENTATION 

Most sandpile models studied so far are characterized 
by irreversible dynamics, where one only adds grains 
which later leave the system through open boundaries. 
Thus, if by adding a grain to a stable state X and relaxing 
it one arrives at another stable state Y, there is usually 
no way to return from y to AT in just one step. As this 



violates the detailed-balance condition and hence drives 
the system out of equilibrium, in our simulations we also 
use the method of Ref. to construct the reversed pro- 
cess that transforms Y back to X. This method is based 
on removing a grain from an (arbitrary) node and then 
fully relaxing the resulting state. 

Following the standard Metropolis algorithm "25, '23l| , 
in each time step we construct a trial state by picking 
at random a lattice node j and adding or removing a 
grain at it. If this renders the configuration unstable, 
an avalanche is generated and the system is fully relaxed 
to a stable state. We then calculate the energy of the 
trial state, E' , compare it with the energy of the origi- 
nal state, E, and accept the trial state with probability 
min(l, ex-p[— {E' — E) / kBT]) . This, however, may lead to 
problems with ergodicity. Since the phase space of our 
model is limited to recurrent states of a Markov process 
defined by the BTW model, there is no doubt that at 
least in theory the system is ergodic. However, some 
transition probabilities can be extremely small. This 
problem can be particularly serious for transitions re- 
lated to large avalanches, as in their case the factor 
exp[—{E' — E)/kBT] can quickly tend to as L — > 00. 

Calculation of a trial state in our model is extremely 
time-consuming, for it requires to perform full relaxation 
of an excited state, and the average number of individual 
topplings in an avalanche is proportional to the number 
of lattice nodes 0- For this reason we were able to 
perform simulations only for rather small lattices with 
N < 9600 (for a pure Ising or BTW models this num- 
ber could be easily increased by a factor 1000). Due to 
the non-local nature of sandpile dynamics employed to 
generate trial states, there seems to be no way to apply 
here any of the many techniques for speeding up Monte- 
Carlo simulations '23'|. The only 'trick' we used was to 
store, for a given configuration, energies of any rejected 
Monte-Carlo trial steps to avoid multiple relaxations of 
the same state. For low temperatures this can speed up 
calculations by a factor of 100. 



IV. RESULTS 

A. Finite-size effects induced by the BTW 
component 

Two essentially different types of finite size effects can 
be expected to appear in the model. The first kind, typ- 
ical of critical points, results from the fact that at criti- 
cality the diverging correlation lengths exceed the system 
size. The second kind is inherent to SOC models with 
open boundaries, as the mean concentration of particles 
near the open boundaries tends to be smaller than that 
in the bulk. For a two-dimensional lattice of linear size L 
this phenomenon brings about a finite-size correction of 
order 0{1/L) This finite-size effect has a noticeable 
impact on the thermodynamics of the hybrid model at 
all temperatures. Its magnitude can be appreciated by 
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FIG. 2: Order parameter m as a function of for lattice 
sizes L = 10, ... , 250. The error bars are less than the line 
width. The dashed line represents a quadratic fit and the 
extrapolated value for i ^ oo is 0.33330(4). 

studying the limit of the temperature going to infinity, 
since in this case our model is governed entirely by its 
SOC component. 

Although the finite-size effects in the BTW model were 
investigated in several papers 24, 25], none of them 
referred to the honeycomb lattice. Figure |21 depicts 
m = (sj), i.e. the order parameter in the standard Ising 
model (magnetization per node) , for several values of the 
system size L. As expected, the dependence of this pa- 
rameter on L is quite strong. Upon approximating m{L) 
as a quadratic in 1/L and extrapolating it to the limit of 
L — > oo, we obtained rrioo ~ 0.33330(4), which suggests 
that ttIqq — 1/3. Note that in our model tti > even at 
infinite temperature and in absence of an external mag- 
netic field. This reflects a property of the BTW sandpile 
model that nodes with 2 grains (corresponding to hj =2, 
or Sj — +1) are more probable than those with 1 grain 
(/ij = 1, or = -1). 

Similarly good quadratic fits were found for the mean 
concentrations c^^\ c^^\ and of nodes with hj = 
0, 1 and 2, respectively. In the limit of i cxd we 
obtained cS} « 0.08334(5), c^i' « 0.29166(2), and c^^ « 
0.62500(3). This suggests that = 2/24, c^^ = 7/24, 

and c^"* — 15/24. Note that the exact values of ci^ 
are known for square |2^ and Bethe |2^ lattices, and 
numerical results are available for hypercubic lattices of 
dimension 2 to 5 [2^ . 

B. Effect of SOC constraints on possible phases of 
the Hamiltonian system 

As we have just seen, the high-temperature properties 
of our model are very unusual: while in standard Hamil- 
tonian systems at high temperatures the magnetization 
disappears (the system is in a disordered, paramagnetic 



phase), our hybrid system exhibits a nonzero, positive m 
as T — *■ oo. Moreover, we will show that to > in the 
limit L oo (of an infinite system), at any temperature 
and for any values of the control parameter J. 

Let Nk denote the total number of lattice nodes with 
k grains, k = 0, 1,2. On the one hand their sum is the 
total number of lattice nodes, 

iV2 + A^i + A^o = (3) 

On the other hand, in any allowed configuration these 
numbers satisfy 

2N2 + Ni>Bl, (4) 

where is the number of bonds on a lattice of size L. 
This relation can be justified using the original version of 
the burning algorithm 0, |1| , which guarantees that the 
El bonds can be mapped onto the N nodes in such a 
way that a node with k grains is an image of at least k 
bonds. 

Equation Q implies c^^^ + c^^^ < 1 and, since in our 

case Bl = 9L^ - 3L, equation igj leads to 2c^^^ -I- c^^^ > 
3/2 - 1/2L. Hence 

TO, = cf -c« = 2(24nc«)-3(cf +c«) > -1/L. 

(5) 

As the right-hand-side of this relation tends to as 
L — > oo, we conclude that there can be no phase with 
a negative value of to. This is not a trivial statement, 
as for suitably chosen values of J and T, to can take on 
any value in the range [0, 1], and for finite lattices there 
even exist allowed configurations with m < 0. An ex- 
ample of a system with to < is a lattice with L = 1 
in which one node has a hight hj = 2 (corresponding 
to the spin Sj — +1) and the remaining 5 nodes have 
hj = 1 (i.e. Sj — —1). Since ((SJ applies to systems at an 
arbitrary temperature, it implies that no paramagnetic- 
ferromagnetic phase transition at finite T is possible in 
the hybrid model, so that the parameter m cannot be 
regarded as a usual magnetization. 

Similar arguments lead to cf' < l/A+l/iL, which im- 
plies that no quadrupolar phase is possible in the model 
(in a quadrupolar phase half of the nodes have spins 
Sj = 0, i.e. c° > 1/2). Consequently, only two phases 
can exist: a high-T paramagnetic (P) and a low-T an- 
tiferromagnetic (A) phase. Note that this reduction of 
the number of possible phases is a consequence of the 
constraints imposed on the Ising model by the sandpile 
dynamics, and has nothing to do with the "criticality" of 
the latter. 

These conclusions are in accord with symmetry con- 
siderations. Because of the constraints imposed by the 
underlying SOC model, the hybrid model has a single, 
unique "ferromagnetic" ground state, and hence the low- 
T system with ferromagnetic coupling (J > 0) is not 
invariant under the Z2 symmetry of the standard Ising 
model; consequently, no ferromagnetic-like phase tran- 
sition is expected in the system. On the other hand. 
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the SOC constraints do not affect the antiferromagnetic 
(J < 0) ground state nor none of its "neighboring" (i.e. 
sHghtly excited) states. Thus, the hybrid model with 
antiferromagnetic couphng at low T obeys the antiferro- 
magnetic Z2 symmetry, which opens a possibility of the 
antiferromagnetic phase transition. 



C. Transition between paramagnetic and 
antiferromagnetic phases 

Since we expect the antiferromagnetic phase to appear 
only for J < 0, we decided to look for a temperature- 
driven transition between the P and A phases only for 
J < 0. For the sake of simplicity, we measure energy 
in units of | J| and set the Boltzmann constant fc^ = 1. 
To investigate a phase transition, we formally divide the 
lattice into two sublattices (such that no adjacent nodes 
belong to the same sublattice) and we define a param- 
eter a as half of the difference between magnetizations 
on each sublattice. In the ground state our system has a 
perfect antiferromagnetic order, with one sublattice com- 
pletely filled with spins s = -1-1, and the other one with 
s = — 1. Hence the mean energy per node u = —1.5, 
magnetization m = and the order parameter a = ±1. 
We expect that \a\ > in the A phase, and a = in the 
high-T phase. 

We start each simulation from the ground state in 
phase A^ with a = 1, and then slowly warm the sys- 
tem up. After equilibrating, at each temperature we use 
an algorithm of [2J| to generate 10* configurations. The 
results are averaged over 10 independent runs. 



1. Internal energy and specific heat 

In Fig. 13 we present the temperature dependence of 
the internal energy per lattice node, m, calculated for a 
system of linear size L = 40. As can be seen, the curve 
has a quite large slope at T « 1.5. This suggests that 
the specific heat 



C = du/dT. 



(6) 



may diverge near T = 1.5, which is the first hint for a 
phase transition. A blow-up of this region is presented 
in the inset, together with results for smaller lattices 
(L = 10 and 20). It shows that the curves obtained 
for different lattice sizes cross each other at Tc ~ 1.475 
and suggests that near this temperature du/dT (i.e. the 
specific heat) significantly depends on the system size. 
Note that this value is different from the critical temper- 
ature of the standard spin-1 Ising model on a honeycomb 
lattice, ri"'"s « 1.2 Hi El . 

The temperature dependence of the specific heat, cal- 
culated for three different lattice sizes from energy fluc- 
tuations, 




FIG. 3; Energy per lattice node, u, as a function of temper- 
ature for J < and lattice size L = 40. The inset presents 
the results for L = 40 (•), L = 20 (□), and L = W (o) in the 
vicinity of T « 1.475. 



is depicted in Fig. 01 As can be seen, C develops a max- 
imum near Tc ~ 1.475, and this maximum grows with 
lattice size L. This L-dependence of C on the system 
size is far more significant near Tc than far from it and 
thus should be attributed to the classical phase transition 
rather than to the effects induced by the SOC constraints. 
The finite-size scaling theory predicts that if there is a 
second-order transition at Tc, C(Tc) should grow as L"/", 
with a, ly being the usual scaling exponents. The depen- 
dence of C at T = 1.475 on L is shown in Figure |S1 as a 
log-log plot. It suggests that the slope tends to 0, which 
implies a = 0. This, in turn, suggests that the hybrid 
system could belong to the same universality class as the 
standard Ising model, for which a = 0, 1^ — 1, and the 
specific heat diverges logarithmically. To verify this pos- 
sibility, in Fig.|31we also plotted a line which represents 
the best-fit of the data to the ansatz C{L) « a -|- 61n(L). 
The fit turns out to be rather convincing. 

Notice the inset in Fig. Q. It depicts C{T) on a log- 
log plot. It can be observed that for large temperatures 
T the specific heat decays as T^^, as it should do for 
usual Hamiltonian lattice spin models. 

A well-established numerical method for analysis of 
critical points consists in studying the temperature de- 
pendence of the fourth cumulant of u, 



Vl = 1- 



(8) 



C 



N 
2^2 



(7) 



This quantity is supposed to have a local minimum at 
Tc, both for continuous and discontinuous phase tran- 
sitions 23, 29]. For first-order transitions the depth of 
this minimum, 2/3 — V^™, carries information about the 
latent heat. The dependence of Vl on temperature in 
our model is presented in Fig. |3 The curves develop a 
clear minimum near T = 1.475. Since its depth decreases 
with L, it suggests that there is no latent heat and the 
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FIG. 7: Specific heat per node, C, for L = 40, calculated from 
energy fluctuations (+ and dotted line) and the derivative 
du I dT ( X ) . The inset depicts the histogram of internal energy 
for T = 10 and L = 40. The solid line is the best-fit normal 
distribution. 



FIG. 5: Finite-size scaling of the specific heat C at T — 1.475 
(log- log plot). The dashed line is the best fit to the formula 
C{L) ^a + bln{L). 



transition is of the second order. 

At this point we apply two tests on u to see if our 
system has some basic properties of standard models of 
statistical physics. First, we verify whether the ther- 
modynamic and microscopic definitions of specific heat, 
Eqs. © and (TJ, are equivalent. The result of this test 
is shown in Fig. |7| which presents C calculated with the 
two methods. The consistency of the results is very good. 
Second, we check if the distribution of energy (far from 
Tc) is normal. An example of a result of this test is pre- 
sented in the inset of Fig. |7| which shows a histogram of 
u as well as the best-fit Gaussian approximation. In this 
case the empirical and theoretical distributions are very 
close to each other. 



2. The order parameter and Binder's cumulant 

The temperature dependence of the antiferromagnetic 
order parameter \a\ is depicted in Fig. |H| It shows that 
the graphs of |a(T)| for various system sizes L cross near 
T — 1.475 (probably at different temperatures). In the 
same region the slope of \a{T) \ gets steeper with increas- 
ing L, while for higher temperatures \a\ quickly tends to 
as L — > cxD. All these properties are typical hallmarks 
of a classical phase transition. However, large error bars 
(not shown) prevented us from doing any reliable quan- 
titative estimates of the critical exponents or Tc based 
on these data. 

We have also performed a study of Binder's cumulant, 
which is the fourth cumulant of the order parameter |2,l| , 

U,^l-^. (9) 
This quantity is particularly useful in finding a critical 
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FIG. 8; The antiferromagnetic order parameter \a\ as a func- 
tion of temperature T for L — 10, 20 and 40. 
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FIG. 9: Binder cumulant Ul as a function of temperature T 
for L = 10, 20, and 40. 



temperature, Tc, which can be extrapolated from the ab- 
scissas of intersection points of Ul{T) for several L, as 
these abscissas are typically only very weakly dependent 
on L. This method is supposed to give precise results 
irrespective of the order of the phase transition . The 
plot of Ul (T) for our model is shown in Fig. As can 
be seen, Binder's cumulants seem to cross each other, 
though the cutting point coordinates are loosely depen- 
dent on L. Moreover, the error bars (not shown), espe- 
cially for larger temperatures, are very high. 

To identify the source of these problems, we investigate 
histograms of the order parameter a at various temper- 
atures. This is motivated by the fact that the theory 
behind Binder cumulant Ul is based on an assumption 
that these histograms can be approximated by a normal 
distribution in a disordered phase, and by a superpo- 
sition of two normal distributions in an ordered phase. 
Four typical histograms are presented in Fig. 1101 The 
first one shows that, for temperature T = 1.2 ^ Tc, the 
histograms are smooth functions of a that can be ap- 
proximated by a Gaussian only for the largest system 
size (L = 40); for small L they develop a slowly decaying 
tail on the left. The second panel shows the situation for 



FIG. 10: The histograms of the antiferromagnetic order pa- 
rameter a for T = 1.2, 1.45, 1.5, 4 and three lattice sizes 
L = 10 (dotted line), 20 (dashed line) and 40 (solid line, nar- 
rowest distribution). 



T = 1.45, which is just below Tc- For L = 10 the system 
can freely switch between states with positive and neg- 
ative a, but such a switch is impossible for L = 40. In 
the intermediate case of L = 20, the system can switch 
to a < 0, but the histogram is very asymmetric, which 
indicates that the actual relaxation time must be much 
larger than the simulation time. As illustrated in the 
third panel, at a little higher temperature T = 1.5 the 
system can already freely switch between a > and a < 
for all investigated system sizes. However, the histogram 
for i = 40 turns out asymmetric about a — Q and very 
erratic, which again indicates that the relaxation time is 
far larger than that available in computer simulations. 
Finally, the last panel shows that at a high tempera- 
ture, T = 4, histograms are again smooth functions of a, 
and each can be approximated by a normal distribution. 
Therefore, for temperatures T « 1.475 and large system 
sizes L, the actual relaxation time is much larger than 
the simulation time. Large relative errors of the data 
presented in Fig. |51 can thus be explained by the phe- 
nomenon of slow relaxation, which seems to be a generic 
property of sandpile models |3Cll |. and which should be 
particularly important at temperatures close to Tc- 



D. Influence of the Ising interactions on criticality 
of the BTW subsystem 

Recall that a signature of criticality is the power-law 
decay for the probability Pki [r) that two nodes at a dis- 
tance r apart have heights k and / (0 < fc, ^ < 3). It was 
proved jlli] that in the bulk of a two-dimensional sandpile 
one has 



-Poo(r) = Pq' + cr-^ + 



(10) 



where Pq and c are some constants that depend on the 
lattice and boundary conditions pTj . This relation can be 
used as a test of the extent to which the self-organized 
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criticality is affected by temperature and interactions. 
Tlie idea is simple: if for T > the system is domi- 
nated by interactions, Pqo (r) should decay in a way typ- 
ical of high-temperature Hamiltonian systems, i.e. expo- 
nentially; if the system is dominated by SOC constraints, 
Poo(r) should decay in accordance with (fTn|l . 

The results for a rather high temperature T = 10, 
where the Ising interactions should not play a significant 
role, are depicted in Fig. ^2 The log- log plot in the inset 
confirms validity of Eq. H10(l . Similar results were also ob- 
tained for all temperatures T > 2.5, whereas for T < 2.5 
an exponential fit was also possible. We explain this as 
follows. Generally one can expect that for temperatures 
T > Tc the decay rate of Pqo (^) will have both power-law 
(SOC) and exponential (thermodynamic) contributions, 
the former proportional to and the latter propor- 

tional to exp(— r/A), where A is the correlation length. 
For r large enough, the power-law should dominate, but 
for smaller distances, like those available in computer 
simulations, Pooif) may be dominated by the exponen- 
tial term, especially for temperatures closer to (but larger 
than) Tc, as in this case the thermodynamic correlation 
length gets large and the exponential term can exceed 
the rapidly decaying algebraic term Another rea- 

son why we did not confirm (|10|l for T < 2.5 is that 
as the system gets cooler, the Ising interactions become 
more pronounced. As one of their effect is to reduce the 
number of sites with hj =0, the statistics for exploring 
-Poo(^) becomes too poor for reliable data analysis. 

Validity of equation (|10|l for T > Tc can be also jus- 
tified by the following physical argument. Suppose that 
the role of interactions, above Tc, is restricted to flipping 
some spins from Sj = to 5^ = ±1, with a probabil- 
ity p, in an uncorrelated manner. This reduces Poo{r) 
by a factor (1 — p)^. Consequently, Poo(?') should sat- 
isfy Hl()|l also for finite temperatures, with the ratio Pq /c 
independent of T. And indeed, our estimates of this ra- 
tio for T = 2.5,3,4,6,8, and 10 yield a constant value 
0.34 ±0.01. 

We conclude that the system preserves the SOC- 
induced correlations at all temperatures T > Tc- Note, 
however, that this has absolutely no impact on the ther- 
modynamic properties of the system. 



V. CONCLUSIONS AND OUTLOOK 

We found that the spin-1 Ising model with nonlocal 
constraints imposed by the Bak-Tang-Wiesenfeld sand- 
pile model, which is expected to exhibit some features of 
both thermodynamic and self-organized criticality, be- 
haves just like a standard model of statistical physics, 
with all hallmarks of a continuous phase transition. Our 
analysis shows that the power-law correlations induced 
by the 'self-organized criticality' of the sandpile model 
have no impact on the order of the phase transition in 
the Hamiltonian subsystem. In particular, we found that 
the critical exponent a « 0, which suggests that the 
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FIG. 11: Probability Pmif) of finding two empty lattice nodes 
r units apart for L = 40 and T — 10. The dotted line is a fit 
calculated from IIOI I. The inset presents the log- log plot of 
|(Po)^ — Poo('')|. The error bars are of order 10~^ (not shown). 



hybrid model belongs to the same universality class as 
the standard Ising model. Similarly, Hamiltonian inter- 
actions of the Ising model cannot destroy "criticality" 
(i.e. power-law correlations) of the BTW sandpile model 
above the thermodynamic Tc- However, unlike thermo- 
dynamic criticality, long-range correlations induced by 
the self-organized criticality do not show up in thermody- 
namic functions (for example, they do not affect temper- 
ature dependence of the specific heat). Clearly, the long- 
range order induced by Hamiltonian interactions has dif- 
ferent impact on the thermodynamics than that induced 
by using the SOC model to reduce the phase space of the 
system. 

Even though the spin and sandpile models we chose 
belong to the simplest models in their kind, their combi- 
nation results in a model that seems intractable mathe- 
matically and is very difficult to treat numerically - hence 
the rather limited system sizes used in the present study. 
It is surely of interest to find a different combination of 
Hamiltonian and SOC models that would exhibit a sim- 
ilar combination of thermodynamic and SOC criticality 
effects, yet would be easier to deal with. 

Compared with the phase space of the standard Ising 
model, the phase space volume of the hybrid model is 
significantly reduced (from 3^ down to with a < 3) 
in a non-local manner. Such reduction resembles that as- 
sociated with lowering of the system dimensionality and 
hence may affect the critical properties of the system. We 
have shown that it actually decreases the number of pos- 
sible phases. Even though our study indicates that the 
critical exponent a in the hybrid model takes on the same 
value as in the Ising model, the question remains whether 
the two systems belong exactly to the same universality 
class and other exponents should be also examined. An- 
other problem is to investigate the decay of correlation 
functions as T | Tc- 

Last, but not least, our numerical results indicate that 
in the BTW model on a honeycomb lattice of linear 
size L — > oo, the probabilities c^'^'' to find k grains in 
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a given lattice site take on particularly simple values, acknowledged. 
= 2/24, cii,' = 7/24, and = 15/24. 
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